Modeling molecules with constraints

When attempting to compute thermodynamic quantities with a molecular simulation, we are frequently confronted with the problem of sampling in a high-dimensional configuration space. The dimensionality of this space is given by the number of degrees of freedom for the molecular system. Techniques which lower the number of degrees of freedom will increase the efficiency of the thermodynamic sampling— provided, of course, that these techniques are physically justified. Thus, an implicit solvent model may be used to eliminate the degrees of freedom associated with the solvent molecules; the standard chemical force fields replace the electron charges with atom-centered partial charges thereby removing the electrons’ degrees of freedom. Further reductions in dimensionality are possible by imposing constraints on the relative positions of the atoms in a molecule. Thus we might specify that the bond lengths and bond angles in a molecule are fixed and only the torsion angles are allowed to vary. It is such a scenario that we examine in this paper. We address two aspects of this problem: how to move a molecule subject to constraints in order to allow equilibrium averages to be computed using the canonical-ensemble Monte Carlo method [1] and how to evaluate the energy efficiently. The imposition of constraints in molecular modeling has been extensively studied [2, §3.3.2, §15.1]. Let us start by elucidating the difference in the treatment of hard constraints in molecular dynamics and Monte Carlo simulations. We treat hard constraints by taking the limit where the “spring constant” for the hard degrees of freedom is infinite. Molecular dynamics simulations then consider the evolution of the resulting system over a finite time. On the other hand, if we wish to determine the equilibrium properties of a system using the Monte Carlo method, then we need to consider averaging over sufficiently long times to allow equipartition of energy among all the degrees of freedom of a system. This is, of